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INTRODUCTION 


Adjerid  and  Flaherty  (refs  1,2)  developed  an  a  posteriori  estimate  of  the 
spatial  discretization  error  in  a  finite  element  method  of  lines  for  solving 
vector  systems  of  parabolic  partial  differential  equations.  They  discretized 
the  system  in  space  using  Galerkin's  method  with  piecewise  linear  finite  element 
approximations.  The  error  estimate  was  calculated  using  Galerkin's  method  with 
piecewise  quadratic  functions.  A  nodal  superconvergence  property  of  the  finite 
element  method  was  used  to  neglect  errors  at  nodes  and,  thus,  improve  com¬ 
putational  efficiency.  Ordinary  differential  equations  (ODE's)  for  the  finite 
element  solution  and  error  estimate  were  then  integrated  in  time  using  the  back¬ 
ward  difference  code  DASSL  (ref  3). 

Adjerid  and  Flaherty  (refs  1,2)  assumed  that  the  temporal  discretization 
error  associated  with  OASSL  was  negligible  compared  to  the  spatial  error.  Thus, 
their  estimate  of  the  spatial  discretization  error  could  be  regarded  as  an  esti¬ 
mate  of  the  total  error.  They  used  their  error  estimate  to  control  mesh  moving 
and  local  mesh  refinement  procedures  that  attempted  to  equidistribute  the  error 
estimate  and  satisfy  a  prescribed  global  error  tolerance.  Similar  mesh  refine¬ 
ment  strategies  have  been  used  by  Bieterman  and  Babuska  (refs  4,5). 

Our  goal  is  to  develop  techniques  that  simultaneously  estimate  the  temporal 
and  spatial  discretization  errors.  To  this  end,  we  consider  M-dimensional  par¬ 
tial  differential  systems  of  the  form 

Ut(x,t)  f(x,t,u,Ux)  =  (D(x,t)Ux(x,t))x  ,  a<x<b,  t>0,  (la) 
subject  to  the  initial  conditions 

u(x,0)  =  uO(x)  ,  a  <  X  <  b  ,  (lb) 


References  are  listed  at  the  end  of  this  report. 


and  boundary  conditions 


AL(t)u(a,t)  +  BL(t)Ux{a,t)  »  gL(t)  , 

(ic) 

AR(t)u(b,t)  +  BR(t)Ux(b,t)  »  gR(t)  ,  t  >  0  . 

The  variables  x  and  t  represent  spatial  and  temporal  coordinates  and  denote  par¬ 
tial  differentiation  when  they  are  used  as  subscripts;  u,  f,  u^,  gL,  and  gR  are 
M-vectors,  and  D,  Al,  Bl,  Ar,  and  Br  are  M  x  M  matrices. 

We,  like  Adjerid  and  Flaherty  (refs  1,2),  discretize  Eq.  (1)  in  space  using 
Galerkin's  method  with  piecewise  linear  finite  elements.  Temporal  discretiza¬ 
tion,  however,  is  performed  by  the  backward  Euler  method  as  opposed  to  using  an 
ODE  code.  A  second  solution  is  calculated  using  trapezoidal  rule  integration  in 
time  and  the  difference  between  the  two  solutions  is  used  to  furnish  an  estimate 
of  the  temporal  discretization  error.  A  third  solution  is  obtained  using 
Adjerid  and  Flaherty's  (refs  1,2)  quadratic  finite  elements  and  the  trapezoidal 
rule  in  time.  This  solution  is  higher  order  in  space  and  time  than  the  original 
piecewise  linear  finite  element-backward  Euler  solution.  Hence,  it  can  be  used 
to  provide  an  estimate  of  the  total  discretization  error  of  the  piecewise  linear 
finite  element-backward  Euler  solution.  Furthermore,  the  difference  between  the 
piecewise  linear  and  quadratic  solutions  calculated  by  the  trapezoidal  rule  can 
be  used  to  furnish  an  estimate  of  the  spatial  discretization  error. 

At  first  sight,  the  above  procedure  seems  to  be  very  expensive;  however, 
nodal  superconvergence  significantly  reduces  computational  complexity.  Defect 
correction  methods  can  also  be  used  to  reduce  costs  associated  with  the  temporal 
integration. 

The  estimates  of  the  temporal,  spatial,  and  total  discretization  errors  of 
the  piecewise  linear  finite  element-backward  Euler  solution  are  used  to  control 
a  global  refinement  procedure  that  attempts  to  keep  an  estimate  of  the  total 
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discretization  error  per  time  step  in  Hi  below  a  prescribed  limit.  Depending  on 


the  proportions  of  the  temporal  and  spatial  error  estimates  to  the  total  error 


estimate,  we  refine  the  time  step,  finite  element  mesh,  or  both. 


The  piecewise  linear  and  quadratic  finite  element  procedures  and  the  tem¬ 


poral  integration  schemes  are  described  in  the  following  paragraphs,  as  well  as 


the  error  estimation  procedures.  Adjerid  and  Flaherty  (ref  6)  proved  that  their 


spatial  error  estimate  converges  to  the  exact  error  as  the  mesh  is  refined  when 


temporal  integration  is  exact  for  linear  parabolic  problems.  Similar  results 


have  not  yet  been  established  when  temporal  errors  are  present;  however,  the 


computational  results  shown  indicate  that  convergence  of  our  temporal,  spatial. 


and  total  error  estimates  are  likely.  Our  global  refinement  strategy  is  then 


presented,  and  it  is  applied  to  an  unstable  heat  conduction  problem.  Finally, 


we  discuss  our  results  and  suggest  some  future  investigations. 


DISCRETE  SYSTEM 


We  simplify  the  presentation  slightly  by  assuming  that  only  Dirichlet  data 


are  prescribed;  thus,  BL(t)  »  BR(t)  =  0,  t  >  0,  in  Eq.  (Ic).  A  weak  form  of 


Eq.  (1)  is  then  constructed  by  multiplying  Eq.  (la)  by  a  test  function 


v(x,t)€Hgl,  integrating  the  result  with  respect  to  x  from  a  to  b,  and 


integrating  the  diffusive  term  by  parts  to  obtain 


(v,ut)  +  (v,f)  +  A(v,u)  =0  ,  t  >  0  ,  for  all  veHg^  • 


The  inner  product  (v,u)  and  strain  energy  A(v,u)  are  defined  as 


b  b 

(v,u)  *  /  vTudx  ,  A(v,u)  »  /  vjDUydx  . 


(2b, c) 


Functions  v  belonging  to  HqI  are  required  to  have  finite  values  of  (v,v)  and 
(Vx,Vx)  and  vanish  at  x  >  a  and  b.  Any  weak  solution  ucHgl  of  Eq.  (2a)  must 


also  satisfy  the  Dirichlet  (essential)  boundary  conditions 


v'' 

75! 
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u(a,t)  »  A[’(t)gL(t)  ,  u(b,t)  «  A-’(t)gR(t),  t  >  0  , 


(2d,e) 


and  initial  conditions  obtained  by  multiplying  Eq.  (lb)  by  v  and  integrating 
with  respect  to  x,  i.e., 

(v,u)  *  (v,uO)  ,  t  s  0  ,  for  all  veHo^  .  (2f) 

A  discrete  version  of  the  weak  system,  Eq.  (2),  is  constructed  by  using 
finite  element-Galerkin  procedures  in  space  and  finite  difference  techniques  in 
time,  as  follows. 

Spatial  Discretization 

In  order  to  discretize  Eq.  (2a)  in  space,  we  introduce  a  partition 

IIr  ;«  |a  =  XQ  <  xj  <  . . .  <  xr  =  b|  (3) 

of  (a,b)  into  N  subintervals  (x^.j^.x-j),  i  =  1,2,...,N,  and  approximate  u  and  v 

by  piecewise  polynomial  functions  U  and  V,  respectively,  with  regard  to  this 
partition.  Thus,  the  spatial ly-discrete  form  of  Eq.  (2)  consists  of  finding 
UcSgNcHEl  such  that 

(V,Ut)  +  (V,f)  +  A(V,U)  =0  ,  t  >  0  ,  for  all  VeSoNcHol  ,  (4a) 

(V,U)  =  (V,u0)  ,  t  *  0  ,  for  all  VeSo^cHol  .  {4b) 

The  spaces  SgN  and  will  be  chosen  to  consist  of  either  piecewise  linear 
or  piecewise  quadratic  polynomial  functions.  The  spaces  of  piecewise  linear 
polynomials  are  denoted  as  and  and  are  easily  constructed  in  terms 

of  the  familiar  "hat"  functions 


4>i(x) 


x  -  x^-.i 


Xi+I  -  X 
i  + 1  "  X  -j 


i_l  <  X  <  Xi 


Xi  <  X  <  Xi+i  ,  i  =  0,1, 


otherwise 


I 


The  piecewise  linear  finite  element  solution  is  written  in  the  form 


.  -  - - .  .  .  .  ^  ^ 
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Ijn  _  (jn-l 

(V.  . -)  +  0C(V.fn)  +  A(V,un)]  + 

Atn 

(l-0)[(V.ff’-l)  +  A(V,Uf’-l)]  =0  ,  for  an  VeSo^  .  (11) 

The  scalar  parameter  0  is  selected  on  [0,1],  ur>(x)  :=  U(x,tn),  etc,,  and  Atp  ;== 
■^n  ■  '^n-1-  simplicity,  the  test  function  V  has  been  assumed  to  be  independ¬ 

ent  of  time,  although  this  will  not  be  strictly  correct  when  refinement  is 
incorporated  into  the  finite  element  method. 

The  two  particular  choices  of  0  that  are  appropriate  for  our  investigation 
are  0=1,  which  yields  the  backward  Euler  method,  and  0  =  >4,  which  yields  the 
trapezoidal  rule.  It  is  well  known  (ref  7)  that  the  local  discretization  error 
of  the  backward  Euler  method  is  0(Atf,2)  a^d  that  of  the  trapezoidal  rule  is 
O(Atn^).  We  will  use  this  difference  in  the  orders  of  accuracy  of  the  two 
methods  to  estimate  the  local  temporal  discretization  error  of  the  finite  ele¬ 
ment  solution. 

ERROR  ESTIMATION 

Local  and  global  estimates  of  the  discretization  error  have  been  success¬ 
fully  used  to  control  refinement  algorithms  that  attempt  to  solve  partial  dif¬ 
ferential  systems  to  prescribed  levels  of  accuracy  (refs  1,2,4-6,8-12).  Our 
goal  is  to  estimate  the  discretization  error  per  time  step  in  solutions  of  Eq. 
(2)  obtained  by  using  piecewise  linear  finite  element  approximations  in  space 
and  the  backward  Euler  method  in  time.  It  seems  most  appropriate  to  gage 
errors 

e  :=  u  -  U  (12) 

in  the  h1  norm 


lleOi  :=  [/  (eTe  +  eTe)dx]’i 
a  ^  ^ 


(13) 
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however,  other  measures  may  also  be  used.  An  error  estimate  that  is  global  in 
space  and  local  in  time  may  at  first  seem  unusual,  but  it  is  commonly  used  when 
spatial  finite  element  approximations  are  combined  with  temporal  finite  dif¬ 
ference  methods  (see  Reference  13). 

Let  the  piecewise  linear  finite  element  solution  obtained  by  using  backward 
Euler  temporal  integration  be  denoted  as  U0El»ri(x)  at  time  tp.  Likewise,  let 
Ut1''^(x)  and  Uj2/r>(x)  denote  solutions  obtained  at  tp,  by  trapezoidal  rule 
integration  with  piecewise  linear  and  quadratic  approximations,  respectively. 

It  is  known  (ref  14)  that  llu(*,tn)  -  Ube^'^(*)IIi  =  ©(Atp^)  +  0(Atn/N). 

Since  llu(*,tf,)  -  L)-[-2,n||^  =  0(Atn3)  +  0(Atp,/N2),  we  should  be  able  to  use  the 
difference  between  Uj2,n  g^d  to  estimate  the  error  in  Ube^''^:  thus, 

Hu  -  Ube^'^i||  4  -  ube^'^Hj  +  iiu  - 

<  -  Ube^'"IIi  ♦  0(Atr,»)  +  0(Atn/N*)  .  (14) 

The  main  problem  in  using  IIUj2»n  -  Ube^'^’IIi  as  an  a  posteriori  estimate  of  Hu  - 
Ube^'^Hi  is  the  computational  effort  required  to  obtain  This  cost  can  be 

reduced  considerably  by  using  the  superconvergence  property  of  the  finite  ele¬ 
ment  method  for  one-dimensional  parabolic  systems.  In  the  present  context, 
superconvergence  implies  that  finite  element  solutions  converge  at  a  faster  rate 
on  than  elsewhere  on  (a,b).  Hence,  the  error  at  the  nodes  may  be  neglected 
relative  to  the  error  in  the  interior  of  the  elements  when  N  is  sufficiently 
large. 

Nodal  superconvergence  has  been  used  by  several  investigators  as  a  means  of 
constructing  a  posteriori  error  estimates  in  finite  element  approximations. 

In  particular,  Adjerid  and  Flaherty  (refs  1,2)  used  it  in  conjunction  with  their 
adaptive  finite  element  method  of  lines.  Their  situation  was  somewhat  more 
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restrictive  than  ours  as  they  also  required  the  temporal  error  to  be  negligible 
relative  to  the  spatial  error. 

The  use  of  the  nodal  superconvergence  property  enables  us  to  approximate 

Uy2,n  as 

(15) 

where  is  obtained  by  solving  Eq.  (7)  using  trapezoidal  rule  integration 

and  ^5  obtained  by  solving  Eq.  (10a)  by  trapezoidal  rule  integration  with 

u2  replaced  by  Eq.  (15).  Furthermore,  it  is  only  necessary  to  test  Eq.  (10a) 
against  functions  v2£§qN,2^  where  is  a  space  of  quadratic  polynomials  that 

vanish  on  Ilfj. 

To  summarize,  our  procedure  for  obtaining  the  finite  element  solution 
UgEl.n  and  its  error  estimate  +  Ej^.n  _  for  the  time  step  [tn-i.t^] 

consists  of: 

(1)  discretizing  Eq.  (7a)  by  the  backward  Euler  method  and  determining 
as  the  solution  of 

Uggl>ri  -  Ube^'^~^ 

(Vl,  . )  +  (Vl,f(*,tn,UBEl'")) 

^^n 

+  A(v1,Ube^'")  =  0  ,  for  all  vleSo^'l  ,  (16a) 

(2)  discretizing  Eq.  (7a)  by  the  trapezoidal  rule  and  determining  Ujl>^  as 
the  solution  of 


UylfH  —  yggl,n'‘l 

(Vl,  . — . )  +  Ji[(Vl,f(*,tn,UTl'"))  +  A(v1,Ut1'") 


At, 


+  (Vl,f(*,tn-i,UBEl'"'M)  +  A(v1,Ube^'"'1)]  *  0  ,  for  all  vUSoN.l,  (16b) 


(3)  discretizing  Eq.  (10a)  by  the  trapezoidal  rule  and  determining  Ej2.n  as 
the  solution  of 
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Uy^*^  +  Ey^*^  “  “  Ey^'^~^ 

(V2.  . -) 

Atn 

♦  5<[(v2,f (..tn.Ujl'"  ♦  Ej^'"))  ♦  A(v2,UTl'"+ET2'n) 

+  (v2,f(..tn-i.UBEl'"‘^+ET2*f'-l))  ♦  A(V2 .Ube^ ' "‘I+Et^ - ^"1 ) ]  =  0  , 
for  all  v2c§qN,2  .  (16c) 

Temporal  error  estimation  is  local;  thus,  we  use  Ube^'*^"^  as  an  initial 
condition  for  the  trapezoidal  rule  integrations  in  Eqs.  (16b)  and  (16c).  Nodal 
superconvergence  and  the  hierarchical  formulation  has  uncoupled  the  piecewise 
linear  and  quadratic  components  of  Uy2,n_  j^e  spatial  error  estimate  £72,0  on 
the  subinterval  (x^.^.x.,-)  is  furthermore  uncoupled  from  the  error  on  other 
subintervals  and  this  significantly  reduces  the  computational  complexity  asso¬ 
ciated  with  solving  Eq.  (16c).  The  solution  of  Eq.  (16t>),  noted  in  step  (2),  is 
necessary  in  order  to  increase  the  temporal  accuracy  of  the  solution  because 
superconvergence  only  increases  the  order  of  accuracy  in  space.  Some  com¬ 
putational  savings  can  generally  be  obtained,  especially  for  nonlinear  problems, 
by  calculating  U-j-l'f’  as  a  defect  correction  to  the  backward  Euler  solution 
Ube^*'’. 

As  described  above, 

en  :>  £72,0  _  Ube^''^*!  (17) 

furnishes  an  estimate  to  the  error  Ru  -  UbeI'*^*!  of  the  backward  Euler  solution. 
Equation  (17)  suggests  the  inequality 

eO  ^  lUjl'O  -  Ube^'^^R^  +  llET2,n||j  _  (10) 

The  term  RUjl'f'  -  Ube^'^"!  ’S  difference  between  two  piecewise  linear  solu¬ 
tions  computed  with  temporal  integration  schemes  of  different  orders  and  can  be 
regarded  as  a  measure  of  the  temporal  discretization  error.  In  a  similar 


manner,  PEy^.niij  can  be  regarded  as  a  measure  of  the  spatial  discretization. 
Indeed,  when  the  finite  element  system,  Eq.  (7),  is  integrated  exactly  in  time, 
Adjerid  and  Flaherty  (ref  6)  proved  that  IIE^tlj  converges  to  the  exact  spatial 
discretization  error  Hu  -  UlPj  as  N  ->  oo  for  linear  parabolic  problems. 

We  conclude  this  section  by  presenting  an  example  that  indicates  that  e^, 
IUjl>n  _  and  llEj^.niij  provide  good  estimates  of  the  total,  temporal, 

and  spatial  discretization  errors,  respectively. 

Example  1:  Consider  the  linear  heat  conduction  problem 

Ut  =  UxxA^  *  0<x<l  ,  t>0,  (19a) 

u(x,0)  =  sin  jrx  ,  0  4  x  <  1  ,  (19b) 

u(0,t)  =  u(l,t)  »  0  ,  t  >  0  .  (19c, d) 

The  exact  solution  of  this  simple  problem  is 

u(x,t)  =  e“t  sin  nx  (20) 

We  solved  Eq.  (19)  on  a  uniform  mesh  with  N  finite  elements  for  one  time 
step  At  using  the  methods  described  above  and  several  choices  of  N  and  At.  The 


effectivity  index 


0  :=  el/llu(»,At)  -  Uggl'lllj 


(see  Reference  15),  is  used  as  a  means  of  gaging  the  accuracy  of  the  error  esti¬ 
mate  el.  Ideally,  we  would  like  0  not  to  differ  appreciably  from  unity  and  to 
approach  unity  as  N  -  «>  and  At  0. 

We  present  a  summary  of  results  for  the  reciprocal  of  the  effectivity  index 
for  a  sequence  of  calculations  performed  with  N  =  ZP+l  and  At  =  5~P/2,  p  = 
0,1,..., 5,  in  Figure  1.  These  results  strongly  suggest  that  0  -  1  as  p  -  <x>. 

We  use  the  temporal  effectivity  index 

0^  •  ,At)  -  Ube^'^'*1 
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as  a  method  of  appraising  the  accuracy  of  the  temporal  error  estimate  - 

Ube^'^“i-  For  fixed  At,  0t  -*  Kt(At)  as  N  -  »  and  the  limiting  value  Kt(At)  -  1 
as  At  -»  0. 

We  solve  Eq.  (19)  for  a  single  time  step  using  a  sequence  of  meshes  with  N 
«  2P,  p  a  3, 4,..., 10,  finite  elements  and  time  steps  of  At  =  0.7,  0.49,  0.343, 
and  0.2401.  We  present  our  findings  for  the  temporal  effectivity  index  0^  as  a 
function  of  p  for  the  four  time  steps  in  Figure  2.  As  expected,  0^  tends  to  a 
limiting  value  'for  large  N,  which  approaches  unity  as  At  -  0. 

Finally,  we  define  the  spatial  effectivity  index  as 

0s  :*  llET2«l|l2/|lu(.,At)  -  Ube^'^Oj  (23) 

and  use  it  as  a  measure  of  the  spatial  error  estimate  IIEj2,l|i2,  For  fixed  N,  0s 
-  K5(N)  as  At  ->  0  and  the  limiting  value  Ks(N)  -►  1  as  N  -►  ». 

Again,  we  solve  Eq.  (19)  and  present  results  for  the  reciprocal  of  the  spa¬ 
tial  effectivity  index  as  a  function  of  At  =  5"P/2,  p  *  1,2,..., 7,  for  meshes 
with  N  a  2,  4,  and  8  finite  elements.  These  results  suggest  that  for  a  fixed  N, 
0s  -  Ks(N)  as  p  -  00,  and  that  Ks(N)  is  reasonably  close  to  unity.  Futhermore, 
it  appears  that  K5(N)  -  1  as  N  increases. 

HESH  REFINEMENT 

The  error  estimates  developed  in  the  preceding  paragraphs  are  used  to 
control  a  simple  global  mesh  refinement  procedure  that  keeps  e^  below  a  spe¬ 
cified  tolerance  TOL.  Suppose  that  a  solution  Ube^'*^"^  and  error  estimates 
Ej2,n-1  and  e^'l  have  been  calculated  at  time  tn-j  using  a  mesh  with  N  elements 
and  time  step  Atp.j.  Further  suppose  that  e*^“l  <  TOL  and  calculate  solutions 
and  error  estimates  at  time  tp  =  tp-i  +  Atp  using  a  mesh  with  N  elements  and 
time  step  Atp  =  Atp-j.  Our  refinement  strategy  consists  of  checking  e'^  and  pro¬ 
ceeding  as  follows: 


(1)  if  ef'  <  TOL,  continue  to  the  next  time  step; 

(2)  if  ef’  >  TOL,  and  0.3  <  IIET2#n|i/en  <  0.7,  double  N,  reduce  Atp  by  30 
percent,  and  redo  the  integration; 

(3)  if  efi  >  TOL  and  0.7  <  double  N  and  redo  the  integration; 

and 

(4)  if  e^  >  TOL  and  llE-r^.nn/gn  ^  o.3,  reduce  Atp  by  30  percent  and  redo  the 
integration. 

Steps  (2)  through  (4)  are  repeated  until  step  (1)  is  satisfied. 

The  main  advantage  of  this  refinement  procedure  is  that  the  separate  esti¬ 
mates  of  the  spatial  and  temporal  errors  allow  different  strategies  to  be  used 
depending  upon  the  dominant  component  of  the  error.  Thus,  if  the  spatial  com¬ 
ponent  of  the  error,  as  measured  by  IIET2'''ll/ef'  is  large,  then  only  spatial 
refinement  is  used  to  reduce  the  total  error.  The  opposite  situation  arises 
when  the  spatial  component  of  the  error  is  small. 

It  is  important  to  note  that  the  error  estimates  used  in  the  refinement 
procedure  are,  at  best,  only  asymptotically  correct.  Thus,  they  will  not  pro¬ 
duce  reliable  estimates  on  coarse  meshes  or  when  errors  are  large.  With  this  in 
mind,  it  may  be  best  to  replace  e^*  by 
refinement  procedure. 

The  specific  choice  of  the  limiting  values  0.3  and  0.7  that  are  used  to 
determine  the  dominant  component  of  the  error  in  our  refinement  procedure  are 
basically  arbitrary.  Under  normal  circumstances,  the  spatial  error  measure 
IIET2,nn/en£(o, 1) ;  thus,  it  is  reasonable  to  divide  (0,1)  approximately  into 
thirds,  i.e.,  (0,0.3),  (0.3, 0.7),  and  (0.7,1)  corresponding,  respectively,  to 
only  temporal  refinement,  temporal  and  spatial  refinement,  and  only  spatial 
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refinement.  This  strategy  may  not  be  appropriate  in  all  situations  and  further 
analysis  and  experimentation  is  needed  to  determine  optimal  refinement  criteria. 

A  local  refinement  strategy,  such  as  those  considered  in  References 
2,4-6,8-12,  is  usually  more  efficient  than  the  global  strategy  presented  herein. 
Our  plans  are  to  combine  refinement  with  a  mesh  moving  method  that  equidistrib- 
utes  a  global  error  measure  on  a  mesh  with  a  fixed  number  of  finite  elements 
(refs  16,17).  It  may  be  possible  to  use  a  simple  global  refinement  strategy  in 
conjunction  with  such  a  mesh  moving  method  since  the  local  error  measure  will  be 
approximately  the  same  on  every  subinterval. 

Doubling  the  number  of  finite  elements  whenever  spatial  refinement  is  per¬ 
formed  simplifies  interpolation  issues,  but  may  add  more  nodes  than  necessary. 
Reducing  the  time  step  by  30  percent  keeps  temporal  accuracy  comparable  to  spa¬ 
tial  accuracy,  since  the  temporal  convergence  rate  is  OCAt^*),  while  the  spatial 
convergence  rate  is  0(1/N).  Thus,  doubling  N  would  correspond  to  reducing  At^, 
by  I//2,  which  is  approximately  30  percent. 

We  apply  the  above  refinement  procedure  to  the  following  singular  parabolic 
problem. 

Example  2:  Consider  the  partial  differential  system 

u^  +  u/[2(l-t)]  =  -Uxx/^  '  0<x<®  ,  0<t<l,  (24a) 

u(x,0)  =  e"X^  ,  0  4  X  <  00  ,  (24b) 

Ux(O.t)  »  lim  Ux(s,t)  =0  ,  0  <  t  <  1  .  (24c, d) 

s-00 

The  exact  solution  of  this  problem  is 


(24c, d) 


u(x,t)  *  .  (25 

This  problem  was  motivated  by  our  interest  in  solving  the  nonlinear 
Schrodinger  equation  in  cylindrical  coordinates  (ref  18).  It  is  known  (ref  19) 
that  the  solution  of  the  Schrodinger  equation  can  "self-focus,"  i.e.,  its 
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solution  can  become  infinite  at  one  point,  while  decaying  elsewhere  on  the 
domain.  Problems  of  this  type  occur  in  laser  optics. 

Such  problems  are  difficult  to  solve  by  traditional  numerical  methods  and 
illustrate  the  need  for  adaptive  strategies.  The  model  given  by  Eq.  (24)  was 
developed  as  a  simple  approximation  of  the  behavior  of  the  Schrddinger  equation. 
Its  solution  "focuses"  in  the  sense  that  u(x,t)  -►  0,  x  >  0,  as  t  -*  1,  while 
u(0,t)  *  1. 

This  problem  was  solved  for  values  of  TOL  of  0.2,  0.1,  and  0.05  and  a  sum¬ 
mary  of  the  results  are  presented  in  Tables  I,  II,  and  III,  respectively.  These 
tables  present  the  relevant  refinement  data  for  each  time  step  of  the  solution 
process.  The  time  and  numerical  parameters  At  and  N  at  the  beginning  of  a  solu¬ 
tion  step  are  found  in  the  columns  labelled  "Initial  Time,"  "Initial  At,"  and 
"Initial  N,"  respectively.  The  refined  values  of  At  and  N,  necessary  to 
complete  the  solution  step  within  the  given  tolerance,  are  found  in  the  columns 
labelled  "Refined  At,"  and  "Refined  N,"  respectively.  The  resulting  time  at  the 
end  of  a  solution  step  is  found  in  the  column  labelled  "Final  Time."  The  last 
column,  labelled  "Total  Error  Estimate,"  lists  the  value  of  ei^  at  the  successful 
completion  of  a  time  step.  The  rows  of  each  table  outline  the  solution  process 
as  it  advances  through  time. 

These  results  indicate  that  it  is  sometimes  possible  to  reduce  the  total 
error  by  refining  only  in  space  or  only  in  time,  and  that  the  error  estimates 
ef’,  -  UbeI'^IIi,  and  can  be  used  to  detect  when  these  situations 

arise. 
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TABLE  I.  NUMERICAL  PARAMETERS  AT  THE  BEGINNING  AND  THE  END  OF 
EACH  TIME  STEP  FOR  SOLVING  EXAMPLE  2  WITH  TOL  =  0.2 


Initial 

Time 

Initial 

At 

Initial 

N 

Refined 

At 

Ref i ned 

N 

Final 

Time 

Total  Error 
Estimate 

0.0000 

0.1250 

16 

0.1250 

16 

0.1250 

0.0885 

0.1250 

0.1250 

16 

0.1250 

16 

0.2500 

0.1496 

0.2500 

0.1250 

16 

0.0875 

32 

0.3375 

0.1062 

0.3375 

0.0875 

32 

0.0875 

32 

0.4250 

0.1931 

0.4250 

0.0875 

32 

0.0613 

64 

0.4863 

0.1381 

0.4863 

0.0613 

64 

0.0429 

128 

0.5291 

0.0782 

0.5291 

0.0429 

128 

0.0429 

128 

0.5720 

0.1144 

0.5720 

0.0429 

128 

0.0429 

128 

0.6149 

0.1404 

0.6149 

0.0429 

128 

0.0429 

128 

0.6578 

0.1911 

0.6578 

0.0429 

128 

0.0300 

256 

0.6878 

0.1052 

TABLE  II.  NUMERICAL  PARAMETERS  AT  THE  BEGINNING  AND  THE  END  OF 


EACH  TIME  STEP  FOR  SOLVING  EXAMPLE  2  WITH  TOL  «  0.1 


Initial 

Time 

Initial 

At 

Initial 

N 

Refined 

At 

Refined 

N 

Final 

Time 

Total  Error 
Estimate 

0.0000 

0.1250 

16 

0.1250 

16 

0.1250 

0.0885 

0.1250 

0.1250 

16 

0.0875 

32 

0.2125 

0.0496 

0.2125 

0.1250 

32 

0.0875 

32 

0.3000 

0.0615 

0.3000 

0.0875 

32 

0.0613 

64 

0.3613 

0.0390 

0.3613 

0.0613 

64 

0.0613 

64 

0.4225 

0.0549 

0.4225 

0.0613 

64 

0.0429 

64 

0.4654 

0.0604 

0.4654 

0.0429 

64 

0.0300 

128 

0.4954 

0.0340 

0.4954 

0.0300 

128 

0.0300 

128 

0.5254 

0.0528 

0.5254 

0.0300 

128 

0 . 0300 

128 

0.5554 

0.0847 

0.5554 

0.0300 

128 

0.0210 

128 

0.5764 

0.0782 

TABLE  III.  NUMERICAL  PARAMETERS  AT  THE  BEGINNING  AND  THE  END  OF 
EACH  TIME  STEP  FOR  SOLVING  EXAMPLE  2  WITH  TOL  =  0.5 


Initial 

Time 

Initial 

At 

Initial 

N 

Refined 

At 

Ref i ned 

N 

Final 

Time 

Total  Error 
Estimate 

0.0000 

0.1250 

16 

0.1250 

32 

0.1250 

0.0451 

0.1250 

0.1250 

32 

0.0875 

64 

0.2125 

0.0261 

0.2125 

0.0875 

64 

0.0875 

64 

0.3000 

0.0377 

0.3000 

0.0875 

64 

0.0613 

64 

0.3613 

0.0375 

0.3613 

0.0613 

64 

0.0429 

128 

0.4041 

0.0211 

0.4041 

0.0429 

128 

0.0429 

128 

0.4470 

0.0300 

0.4470 

0.0429 

128 

0.0300 

256 

0.4470 

0.0176 

0.5070 

0.0300 

256 

0.0210 

256 

0.5280 

0.0202 

0.5280 

0.0210 

256 

0.0210 

256 

0.5490 

0.0272 

DISCUSSION 

We  developed  methods  for  calculating  a  posteriori  estimates  of  the  total, 
spatial,  and  temporal  discretization  errors  when  a  vector  system  of  parabolic 
partial  differential  equations  is  solved  using  piecewise  linear  finite  elements 
in  space  and  the  backward  Euler  method  in  time.  The  error  estimates  are 
obtained  by  using  higher  order  methods,  with  nodal  superconvergence  being  used 
to  improve  computational  efficiency. 

The  three  estimates  were  used  to  control  a  global  refinement  procedure  that 
attempts  to  keep  a  global  measure  of  the  error  per  time  step  below  a  prescribed 
tolerance.  Refinement  can  be  performed  in  space,  time,  or  both  space  and  time 
depending  on  the  dominant  component  of  the  error  estimate. 

Comparison  of  the  exact  and  estimated  errors,  presented  in  Example  1,  give 
us  some  confidence  in  the  accuracy  of  our  error  estimates.  Additionally,  the 
results  of  Example  2  provide  an  indication  of  the  utility  of  these  estimates  in 


an  adaptive  procedure.  In  certain  situations,  only  spatial  or  temporal  refine¬ 
ment  was  needed  to  keep  the  total  error  within  the  prescribed  tolerance,  and 
our  error  estimates  could  be  used  to  determine  when  these  situations  arise. 

This  is  the  first  attempt  that  we  know  of  which  simultaneously  addresses 
spatial  and  temporal  errors  with  different  refinement  strategies.  Some 
researchers  (refs  8-11)  have  used  binary  refinement  in  space  and  time,  but  did 
not  attempt  to  determine  the  dominant  component  of  the  total  discretization 
error.  As  noted,  method  of  lines  techniques  (refs  1,2, 4, 5)  typically  assume 
that  temporal  integration  is  exact  and  refine  based  on  estimates  of  spatial 
errors.  There  is  a  great  potential  for  techniques  that  utilize  different  spa¬ 
tial  and  temporal  refinement  strategies,  particularly  with  problems  having 
singularities.  Our  work,  however,  is  still  very  preliminary  and  there  is  still 
a  great  deal  to  be  done.  Rigorous  convergence  results  for  our  error  estimates 
are  yet  to  be  established.  The  refinement  algorithm  previously  discussed  is 
very  simple  and  will  likely  benefit  from  further  experimental  and  theoretical 
analyses.  We  also  anticipate  that  the  inclusion  of  a  mesh  moving  procedure 
based  on  equidistributing  a  global  error  measure  (ref  16)  will  dramatically 
improve  the  performance  of  our  adaptive  solution  technique.  In  the  future,  we 
would  like  to  extend  our  techniques  to  multi-dimensional  problems  and  to  con¬ 
sider  higher  order  spatial  and  temporal  discretization  methods. 
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ciprocal  of  the  spatial  effectivity  index  9^  for  N  =  2 
:urve  1),  N  *  4  (curve  2),  and  N  =  8  (curve  3)  versus  p 
lere  At  =  5-P/2  when  solving  Example  1.  Note  that  Og 
iproaches  1  as  N  increases. 
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